Two-body problem in graphene 



o 

(N 
C 

a 
oo 



J. Sabio 1 ' 2 , F. Sols 2 , and F. Guinea 1 

1 Instituto de Ciencia de Materiales de Madrid ( CSIC), 
Sor Juana Ines de la Cruz 3, E-28049 Madrid, Spain 
Departamento de Ftsica de Materiales, Universidad Complutense de Madrid, 28040 Madrid, Spain. 

We study the problem of two Dirac particles interacting through non-relativistic potentials and 
confined to a two-dimensional sheet, which is the relevant case for graphene layers. The two-body 
problem cannot be mapped into that of a single particle, due to the non-trivial coupling between 
the center-of-mass and the relative coordinates, even in the presence of central potentials. We focus 
on the case of zero total momentum, which is equivalent to that of a single particle in a Sutherland 
lattice. We show that zero-energy states induce striking new features such as discontinuities in the 
relative wave function, for particles interacting through a step potential, and a concentration of 
relative density near the classical turning point, if particles interact via a Coulomb potential. In 
the latter case we also find that the two-body system becomes unstable above a critical coupling. 
These phenomena may have bearing on the nature of strong coupling phases in graphene. 
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I. INTRODUCTION 

The problem of interactions in graphene layers was a 
subject of research even before its isolation and char- 
acterization in the laboratory 1 ' 2 . It is peculiar due to 
its two-dimensional nature and to the honeycomb lattice 
structure into which ions are arranged. In the low-energy 
limit, the electronic properties are described by a Dirac- 
like equation for massless and chiral electrons^. Undoped 
graphene has a vanishing density of states at the Fermi 
level and therefore a diverging screening length, so in- 
teractions are expected to yield a more singular behav- 
ior than in a conventional Fermi liquid picture, already 
well established in metals and electron gases^. In fact, 
a weak-coupling scaling analysis was early performed^ 7 -, 
shedding light on the role of electron-electron interac- 
tions mediated by the Coulomb potential: the resulting 
divergences in perturbation theory, once handled conve- 
niently, turn out to have a small effect on the low energy 
properties of electrons, as they are marginally irrelevant 
in the renormalization group sense. However, that analy- 
sis does not exclude the possibility of phases with broken 
symmetry, where interactions play a major role, when 
the dimensionless coupling g = e 2 /ehvp is of order unity 
or larger, e being the electronic charge, vf the Fermi ve- 
locity, and e the dielectric constant of the environment in 
which graphene is embedded 8 - - —. This is expected to be 
the relevant case for graphene layers in vacuum, where 
e = 1 and g ~ 2.16. Along that direction, several works 
in the recent literature are pointing out the existence of 
exotic strong coupling phases in graphene layers, where 
pairing of electron and holes could give rise to an exci- 
tonic state, where a gap would be opened rendering the 
system insulating 8 - ' 11 ' 12 . 

On the other hand, it is common wisdom in the field 
of strongly correlated systems, that the study of the in- 
teraction between two particles can provide important 
insights on the many-body physics. The relevance of this 
kind of studies in graphene has already been shown when 



addressing the Coulomb impurity problemi^r— . Here, a 
critical value of the coupling marks the breakdown of the 
Dirac vacuum, whose study requires consideration of the 
whole many-body problem. As pointed out recentl y 17 ' 18 , 
there could be a relation between this instability and 
the formation of an excitonic condensate in the strongly- 
coupled many-body problem. Thus two-particle physics 
seems to underlie many features of the full many-body 
problem. 

In this paper, we address the problem of two interact- 
ing Dirac electrons in two spatial dimensions, mediated 
by non-relativistic central potentials. This feature makes 
the problem different from the already well addressed 
fully relativistic problem. One of our goals is to shed 
light on the relation between the two-body problem and 
the many-body instabilities, so we pay special attention 
to the case of the bare Coulomb potential. However, 
the general problem happens to show peculiarities that 
make its separate study worthwhile. Importantly, the 
two-body problem cannot be mapped exactly into the 
one-body Coulomb impurity problem on the same lattice. 
In fact, we will show that the two-body problem presents 
remarkable differences, the most important being the sin- 
gular role played by localized zero-energy states at those 
points where the kinetic energy vanishes. Interestingly, 
we find that, for zero center-of-mass momentum, the two- 
body problem in graphene is equivalent to that of a single 
particle in the lattice model proposed by Sutherland 19 , 
where both lattices are considered in their continuum 
limit. The presence of zero-energy states can induce non- 
analyticities in the relative wave function, giving rise to 
partial localization phenomena for the Coulomb interact- 
ing case. In order to get a better understanding of the 
novel features, we address first the case of two particles 
interacting via a step potential. 

The paper is organized as follows: Section II presents 
some general features of the two-body problem. Section 
III addresses the case of zero center-of-mass momentum, 
which is simpler to analyze and of potential relevance to 
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many-body instabilities. Sections IV and V study the 
case of step and Coulomb potentials, respectively. Sec- 
tion VI discusses some aspects of the case of arbitrary 
center-of-mass momentum. Section VII is devoted to a 
discussion of the relation of the present work to many- 
body phenomena. Finally, section VIII summarizes the 
main conclusions. The paper includes three appendices 
where some technical issues have been collected. 



II. GENERAL FEATURES 

As the main applications of this problem concern 
graphene sheets, we start with a formulation of the prob- 
lem in terms of the continuum theory of graphene elec- 
tron motion, which is known to be described by the Dirac 
equation. For a single particle, the wave function is a two- 
component spinor characterized by the quantum numbers 
of spin and valley, both with degeneracy two. For zero 



magnetic field and scalar translationally invariant inter- 
actions such as those which we will consider here, we 
can neglect the effect of those extra degrees of freedom 
and concentrate on the two-component spinor descrip- 
tion. Later, we will discuss the effect of the spin. The 
single particle Dirac equation reads: 



v(r) 
-id x + d y 



-id x - d y 
v(r) 



*A(r) 



= E 



*a(f) 



(1) 



where the pseudo-spin index i = A,B refers to the two 
inequivalent sites within the unit cell of the honeycomb 
lattice. 

Now let us consider the two particle problem. We con- 
struct two-particle wave functions from the tensor prod- 
uct of single-particle ones, ^(ri, r 2 ) = ^i(ri) ® \E , J (r 2 ). 
This allows us to write the Schrodinger equation for 
the interacting problem, that in the language of four- 
component spinors reads: 
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As we are dealing with translationally invariant poten- 
tials, we can switch to the center-of-mass frame, defining 
the new coordinates R = ri j, 1 " 2 and r = ri — r 2 . It is 
also convenient to apply the unitary transformation 



and use a plane wave ansatz for the center-of-mass part 
of the wave function, ^j(R,r) = e lK ' R ^j(r). We arrive 
at the following eigenvalue problem: 
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where 9k = asctaxi{K y / K x ) and polar coordinates are 
used for the relative coordinate. As a first remark on 
this equation, we notice that the center-of-mass coordi- 
nate does not decouple from the relative one, even though 
the potential only depends on the latter. This is a conse- 



quence of the chiral nature of the electron carriers, where 
pseudo-spin and momentum are coupled. This kind of 
coupling also prevents the Hamiltonian from commuting 
with the relative angular momentum, thus frustrating a 
possible decomposition of the problem in terms of partial 



III. THE CASE K = 



In order to gain insight into the many-body problem, 
the most interesting case is that of zero total center-of- 
mass momentum. Then the two particles have opposite 
momenta, like in the Cooper channel in metals. Any 
pairing effect should be particularly important in this 
energetically most favorable case. It is also the simplest 
one, because it decouples the second component ^> 2 (r) 
from the rest. In effect, the equation for this component 
reads 



[v(r) -E]4> 2 {r) = 



(5) 



whose solution is ip2{r) = except at the particular 
point v(r) = E, if it exists. That point having measure 
zero, we will ignore the ip2 component as physically irrel- 
evant. However, we will see later that, at the point where 
v(r) = E is satisfied, zero-energy states are responsible 
for important non-analyticities in the other components. 

We are thus left with an effective three-component 
problem. Remarkably, the K = Hamiltonian commutes 
with the relative angular momentum, so we can use the 
following ansatz for the wave function: 



s<('+i)*03 (r) 



(6) 



where the prefactors have been chosen for convenience. 
The labeling of the components in Eq. (JSJ) has been 
changed in order to accommodate it to the three- 
component case. After using this ansatz, the system of 
equations reads: 



FIG. 1: Scheme of the lattice proposed by Sutherland in Ref. 

The two-body problem in the low-energy sector of the 
honeycomb lattice, for K — 0, is mathematically equivalent 
to a single-particle problem in this lattice. Zero-energy states 
in the Sutherland lattice appear due to the existence of a flat 
band. 



A. Symmetry properties 

Let us now analyze the symmetry properties of the 
K = solutions. In the original basis, the two-body 
wave functions read: 
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*BA(ri,r 2 ) 
* BB (ri,r 2 ) 
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i(r) 



e 
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where the symmetric combination has been taken i/; 2 = 0, 
as argued above. This wave function has a spinorial 
structure, due to the pseudo-spin of the particles, and 
a spatial structure coupled to the former. The symme- 
try properties under exchange of particles are studied by 
doing the transformation: 



v(r) -E B r + l - 
-2(9, - v(r) - E 2(d r + *±i) 
-Br + l - v{r) - E 



Mr) 
Mr) 
Mr) 



(7) 

It is interesting to note that these equations, as well 
as those directly derived from the full Hamiltonian for 
K = 0, Eq. (f4]), can also be obtained as the contin- 
uum limit of a one-particle lattice Hamiltonian, defined 
in a triangular lattice with three sites in the unit cell, 
as initially considered by Sutherland 1 —. A scheme of this 
lattice is shown in Fig. [TJ 

Within this formulation, the case I = is the most 
symmetric one: 



e-^Mr) 
e^Mr) 



(8) 



Henceforth, we will refer to it as the s-wave, and as we 
will see, simple solutions can be obtained for this case 
taking advantage of its symmetry. 



ri ^ r 2 



(10) 



The first transformation, for K — 0, translates into 
cj> (j> + 7T. It follows inmediatly that wave functions 
with / even are antisymmetric under particle exchange, 
while those with I odd are symmetric. Hence, the s-wave 
is, interestingly, antisymmetric. This somewhat counter- 
intuitive result reflects the role of the sublattice pseudo- 
spin in the orbital wave function. 

This has consequences on the total wave functions, 
once both spin and valley degrees of freedom are also 
considered. Everywhere in this article we assume that 
the two particles belong to the same valley. Since their 
total wave function must be antisymmetric, the following 
two families of solutions appear: 



(i) *jc=o,i=odd(ri,r 2 ) ® ^ (| U) ~ I It)) 

(ii) **r=o,i=even(ri,r 2 ) ® | ff) 

*K=0,i=even(ri,r 2 ) (g) j Lj.) 

*jf=0 l I= O v TO (ri,r 2 ) <8> ^ (| U> + I ID) 



(11) 



4 



Therefore, the lowest angular-momentum channel (I = 0) 
corresponds to a triplet spin-state, as opposed to what 
happens with ordinary Schrodinger electrons. 



B. Other mathematical properties 

Equation (J7J comprises three coupled differential equa- 
tions. Adding the first and the third equations we may 
solve for fa in terms of fa and fa: 









k,= 


!V E)/2 



E/2 



02 = ^7 £ ( r )(01 + 0s) , 



(12) 



where e(r) = E — v(r). Substracting the same two equa- 
tions, we obtain 



d r fa = e(r)(fa - fa) . (13) 
T2|) and relating the result to Eq. 



Differentiating Eq 
(fill)) , we obtain 



dr[e{fa + fa)} = - [{21 - l)fa - {21 + 1) 
r 



(14) 



On the other hand, the second equation of system ([7]) 
can be rewritten as: 



OA 







I - 1 



4/ 



1 + 1 



4/ 



(15) 

Thus, system (JTJ) can be solved in principle by first solv- 
ing for fa and fa from (fT4|) and ([To]) and then obtaining 
fa from (fl"2j) . The absolute values of fa and fa should 
remain bounded as long as e{r) is bounded. A problem 
may appear for r — ¥ in the case of the Coulomb poten- 
tial (v(r) ~ r _1 ). In such a case, we will see that regular 
solutions can be obtained analytically for r — > 0, yielding 
in fact a useful starting point to initiate the numerical 
integration. 

Another important issue arises when e — > 0, i.e., at 
those points where the kinetic energy vanishes, when- 
ever I > 0. For a smooth potential, a linear approx- 
imation of e(r) around the vanishing point ro holds, 
e{r) = Xx + 0{x 2 ), where x = r — ro. The differen- 
tial equations (|14[) and (|15p can be approximated around 
this point, yielding: 



d ( 

rfa; 



03) 



+ ^,)] ~ 
1 

~ ( 



+ 03) 



(16) 
(17) 



FIG. 2: Scattering of two particles interacting through a short 
range potential. 



to rp in fa and fa, while fa will remain continuous: 



CH log(ar) 

r a; 



C 2 



2/ 

-Ci log (x) 

ro x 



(18) 



Notice, however, that these non-analyticites give a finite 
contribution to the probability J drr \4>\ 2 - Therefore they 
are physical solutions of the Dirac equation. 

We end by noting that this singular behavior remains 
unaltered even in the pres- ence of a small mass in the 
two-electron Dirac equation. Mathematically, this hap- 
pens because the mass terms in the equations are sub- 
leading in the short-distance expansion around the point 



IV. STEP POTENTIAL 

Some physical insight into the subtle properties of the 
interacting two-particle problem can be obtained by con- 
sidering the simpler case of a step potential, which is 
typically considered a good effective description of the 
more general class of short-range potentials: 



v(r) = 



v r <r 
r > ro 



(19) 



As usual, the procedure is to construct the solutions 
for each region and eventually match them, as shown 
schematically in Fig. [2] 

For arbitrary energy E, the solutions are given by 
Bessel functions of the form: 



The solution for these equations reads fa + fa ~ —2Cijx 
and fa — fa ~ 2Ci + (2C 2 /r ) log(x). We see that, for 
I > 0, a smooth potential will show non-analyticities close 



Xi 



aJi-i(kr) 
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where the coefficients and the eigenvalues are determined 
from the diagonalization of Eq. ([7]), the result being 

E = vo + 2k, (i, 1, i ) (21) 
E = v - 2k, ( i, -1, i ) (22) 

£ = wo, ^= ( 1, 0, -1 ) (23) 

the normalization being chosen to simplify the global 
wave function, Eq. (|6|). The first solution corresponds 
to two electrons located in the upper Dirac cone, while 
in the second solution the two electrons are in the lower 
cone. The third solution describes the case where one 
particle is in the upper cone and the other one in the 
lower cone. Due to the zero total center-of-mass momen- 
tum, this solution has zero total energy. Notice that, for 
fixed E, the relation of k with the energy depends on the 
solution chosen. The third one is valid for arbitrary k. 
Importantly, when E — vq there are other zero-energy 
states that are also solutions of the two-particle Dirac 
equation. They have the form 

X3 = ^ | ( a-Z+l )2]l/2 ( ' a+l+1 ) ' ( 24 ) 

where a is a continuous parameter that can take any real 
value. These polynomial solutions are in general non- 
physical, as they cannot be properly normalized. How- 
ever, they can be considered responsible for the non- 
analyticities shown to exist at those points where the ki- 
netic energy vanishes (see section III.B). For the step po- 
tential, which is non-analytic itself, the existence of zero- 
energy states induces discontinuities in the radial wave 
function, thus changing the usual matching conditions. 
How such an anomalous behavior arises is explained in 
detail in Appendix A. Notice that in the presence of a 
mass this kind of solutions still exist, since they involve 
large derivatives within a narrow distance range. 

We note that the traditional criterion of imposing con- 
tinuity of the wave functions does not work here due to 
the insufficient number of matching parameters. To ex- 
emplify this issue, let us consider the situation shown in 
Fig. [2J where ki = (v Q - E)/2 > and k u = E/2 > 0. 
The solution for region I only includes Bessel functions 
of the first kind, J;(fcjr), as those of the second kind 
ones are singular at the origin. Hence <fi( ~ AiJi^kir). 
For region II both solutions must be considered, (pj 1 ~ 
BiJh (kir) + BiYii (kjj(r). Due to the freedom for global 
normalization, only two relative values of the three con- 
stants are relevant. Thus only two parameters remain 
to satisfy the continuity of three equations, one for ev- 
ery component of the wave function, leaving the problem 
overdetermined. 

In Appendix A it is shown that, when zero-energy 
states at the matching point are taken into account, a 
third matching parameter arises naturally which permits 
a discontinuity in the radial wave function. Namely, we 
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FIG. 3: Numerical solution of the step-potential for K = 
and I = 1 (roE — 10, E = Vo/2). Top: first component of 
the radial wave function. The discontinuity induced by the 
zero energy states arises naturally in the numerical solution. 
Bottom: second component of the radial wave function. As 
also predicted by the new matching conditions, the second 
component does not have a d iscontinuity 



obtain 

A^M^VoWiM = -2C 

A0 2 (r o ) = 

A0 3 M = -2C, (25) 

where C is the extra parameter to determine. With the 
new matching conditions, an exact solution of the K = 
case becomes possible, as detailed in Appendix B. It is 
also interesting to compute the solution of the differential 
equations numerically, as shown in Fig. [3] Here, the first 
two components of the radial wave function are plotted. 
As predicted in equation (I25|) , the first component shows 
a discontinuity induced by zero-energy states at the point 
r = tq. The same consideration applies to the third 
component (not shown) but not for the second one. 

It is also worth noting that, as expected, the s-wave 
shows a simpler behavior by virtue of its symmetric form. 
In this case, inspection of Eq. ||7J) shows </>i(r) = — fefr) 
and the problem reduces to a two-component one. The 
matching conditions reduce then to continuity and the 
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s-wave problem essentially behaves as that of the single 
particle. As a corollary, Eq. (f25j) leads to C = in this 
case. Thus we may state that the I = case has a struc- 
ture similar to that of the impurity one-body problem in 
graphene. 

V. COULOMB POTENTIAL 

We now turn to the more relevant case of a long-range 
Coulomb potential, v(r) = g/r, where g = e 2 /ehvF for 
low-energy graphene electrons. It is convenient to find a 
more suitable form of Eq. ([7]) in order to obtain analyt- 
ical solutions when possible. This is done by the usual 



procedure of analyzing the short and long distance limits. 
At short distances the wave function components have 
the form <fo(r -> 0) ~ r'' 1 / 2 , with 7 2 = \(l + Al 2 - g 2 ). 
On the other hand, the long-distance wave function be- 
haves like a plane wave of the form <f>i(r — > oo) ~ e ±lB ' r ' 2 . 
Hence, we can make the following ansatz for the radial 
wave function: 

Up) = P^e-^Up) (26) 

where we have introduced the dimensionless radial com- 
plex coordinate p = iEr. By applying this transforma- 
tion, Eq. (JTJ) becomes 



ip + g pd p - £ + 7 + 1 - 

-2(pd p -% + 7 -J + I) ip + g 

-pd p + § - 7 + H 



2 

2(pd p -% +7 + ^ + 1) 
■3 *P + .9 



r 



^i(p) 
* 2 (p) 
k(p) 



= o. 



(27) 



The general case is difficult to handle and only the s-wave 
channel admits an analytical solution, since it reduces to 
an effective single particle problem. The details of this 
solution are summarized in Appendix C. For general an- 
gular momentum I, Eq. (|27[) must be solved numerically. 



Before addressing the full solution, let us point out the 
remarkable behavior of the wave functions at short dis- 
tances. As we have seen, it goes like a power law r 7-1 / 2 , 

where 7 = |-\A+ 4 ^ 2 _ 9 2 - For \9\ < 9c = VT+4F 
only 7 > is acceptable. By contrast, for \g\ > g c , the 
7 parameter becomes imaginary and the wave function 
shows a pathological short-distance behavior, going like 
r~ 1//2 [cos(|7| log r) ±« sin(| 7 | logr)]. Thus the wave func- 
tion oscillates dramatically towards the center. This kind 
of behavior was already found in the Coulomb impurity 
problem, where it was related to an instability of the wave 
function that could signal the breakdown of the Dirac 
vacuum. For strong enough couplings, the two-particle 
interaction would produce electron-hole pairs from the 
vacuum, and a full quantum field-theoretical treatment 
of the problem could be necessary. The consequences 
for the two-body problem of this effect have not been 
addressed in this paper, although from the study of the 
Coulomb impurity problem, we may expect a non-linear 
screening as a result of the reorganization of the many- 
body vacuum.— 

In order to gain further insight into the instability, a 
semiclassical analysis like that performed in Ref. LLJ can 
be applied. Our starting point is Eq. ([7]) for the Coulomb 
potential with the ansatz <fn (r) = e ^_ , where p r is a 
slowly varying function of r. This justifies the assump- 
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FIG. 4: Sketch of the Coulomb potential (attractive, in this 
case), including all the relevant lenght scales discussed in the 
text. E is the total energy of the particle, ro is the classi- 
cal turning point, and the region between n and is the 
classically forbidden region. 



tion d r p r — 0. Equation ([7J then reads 



\-E i2 Pr + i=iZ2 
-Upr + ^ *-E 
—i2p r 





i4 Pr + 2l±l 
'+1/2 a. -e 



= 0. 



(28) 



The determinant vanishes if one of the following relations 
is satisfied: 



E. 



5 2 = (--£) 2 -4(4i 2 + l). 



(29) 
(30) 
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FIG. 5: Numerical solution of the radial wave function for the 
case of Coulomb interaction and zero center-of-mass momen- 
tum. The chosen angular momentum is I = 1. Top: wave 
functions for g — —0.5 and E > 0. Bottom: wave functions 
for g = 0.5 and E > 0. Notice that, where the condition 
g/ro = E is satisfied, a singularity is induced by the localized 
zero-energy states. 



The first equation defines a point, r$ = g/E, where any 
function p r gives a solution thanks to the existence of 
zero-energy states, as discussed in the previous section 
for the step potential. Notice that this condition only 
is fulfilled for repulsive electrons with positive energy 
or attractive electrons with negative energy, two cases 
which are related through a symmetry transformation. 
The second equation, on the other hand, defines a non- 
classical region where p% < 0. The region is r\ < r < r2, 
where r\ ; 2 — -§ T ■gv / 4P + 1. Remarkably, we see that 
r\ < tq < rii i.e. the point where zero-energy states 
nucleate belongs to this classically forbidden region, as 
sketched in Fig. |4j 

The Coulomb problem still presents other peculiarities. 
In Fig. [5j a numerical estimate of the radial wave func- 
tions is shown for I = 1 and two different signs of the 
interaction below the critical value. The main feature in 
this solution concerns again zero-energy states: when the 
condition g/r^ — E is fulfilled, zero-energy states must 
be taken into account and become responsible for sin- 
gularities in the first and third components of the wave 
function when I ^ 0, as shown above on general grounds. 




1 2 3 4 5 

rE 

FIG. 6: Probability of finding one electron within a distance r 
from the other electron, P(r) = £V J Q r d 2 r|</>i(r)| 2 , for K = 0, 
1 = 1, and various values of the dimensionless Coulomb cou- 
pling constant g, in the important case where a classical turn- 
ing point exists (rE = g). Here, the influence of zero-energy 
states translates into (i) a suppression of the electron density 
for r < ro, and (ii) a tendency to concentrate the probability 
near r = ro as the critical coupling g c is approached. The 
results are normalized to their value at rE = 5. 



In the Coulomb case, this effect has remarkable conse- 
quences, such as a drastic suppression of the probability 
of finding the particle in r < ro, and a tendency to in- 
creasingly localize the radial wave function near r = ro 
when the critical point g c is approached from below (see 
Fig. [6]). We find numerically that for g = g c the relative 
wave function becomes effectively localized near r = tq. 

These results are similar to those obtained for a sin- 
gle particle in the Sutherland lattice^ with a Coulomb 
potential. As shown in Fig. [7J we use a 30 x 30 lat- 
tice with the structure of Ref. [H (see Fig. [1]) and 
periodic boundary conditions. The potential is v(r) = 
v e~ r / rd I \/r 2 + r\, with w = t > 0,rd — 20a and 
7*i = 0.5a, where t is the hopping, and a is the distance 
between nearest neighbor equivalent atoms. The states 
considered to construct the density plots are in the range 
of energies 0.25i < E < 0.35t. Since they are not eigen- 
states of the Hamiltonian, they are expected to contain 
several angular channels. However, as shown in Fig. [7J 
this energy spread is sufficient to find an enhancement 
of the density in the region near E = v(r). When states 
in the range E < are considered, the density shows a 
delocalized distribution. In the plots, notice that there 
are details coming from the underlying lattice structure 
that are not relevant for our discussion, since we focus 
on the continuum limit. 

In case of considering a small mass in the prob- lem, 
the two most prominent features of the Coulomb two- 
body problem, i.e., the instability above a critical cou- 
pling and the influence of zero-energy states, are not es- 
sentially altered. The Coulomb instability has its origin 
in the short-distances behaviour of the wave function, 
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FIG. 7: Logarithmic density distribution of the two-particle 
wave function interacting with a Coulomb potential (see text 
for details), for K = 0. One of the particles is assumed to 
be placed at the center of the square. The equations are 
discretized using the Sutherland lattice, and the energy range 
of the integrated states is chosen such that they cover the 
region where the condition g = tqE is fulfilled. As expected, 
the results show a clear concentration of the density at this 
point ro. The short-length features of the density plot reflect 
the underlying lattice structure. 



where mass terms are subleading. These terms are also 
subleading in the expasion around the classical turning 
point rO, thus not helping to prevent the appearance of 
non-analyticities caused by zero-energy states which also 
exist for nonzero mass. 



energy states reading now 



r a sm a (6 K - 



-r a sm a {0 K 



-A* 








(31) 



where A\ and A 2 , as well as a and j3, can take arbitrary 
values. A similar analysis to that of K = can be per- 
formed here. For the case of a step potential, they trans- 
late into a change in the matching conditions, since the 
introduction of two new parameters (B\ and B 2 ) changes 
the continuity of the wave function: 



A^i(r o ,0) 
A^ 4 (r 0) ^) 



BxifrK, e K ) 
B 2 (<t>,K,9 K ) 




= -e™«By(^K,B K ) 



(32) 



Linearization of the equations close to the point where 
e(r) = (the classical turning point) shows that, even 
for smooth potentials, singularities in the relative wave 
functions arise. Still, they are square integrable and give 
a finite contribution to the probability. 

As for the Coulomb instability, its existence can be 
probed by checking the short-distance limit of the full 
Hamiltonian given in (j4]), for the case v(r) = g/r. It is 
not difficult to see that the small r limit is controlled by 
if-independent terms. Thus, for r — > we recover the 
K = limit, where an instability has already been iden- 
tified. Hence, we expect this instability to be a general 
feature of the Coulomb problem. 



VII. RELEVANCE TO MANY-BODY 
PHENOMENA 



VI. EXTENSION TO FINITE K 

The most salient features we have found in the problem 
of two interacting particles are, so far, the influence of 
zero-energy states and the appearance of instabilities for 
the Coulomb potential. Next we ponder to what extent 
those results still apply for the general case of nonzero 
center-of-mass momentum. 

Zero-energy states are investigated by taking v(r) = E 
in the eigenvalue problem (j4]) . Inspection of the resulting 
Hamiltonian reveals that it separates into two indepen- 
dent sectors. Hence, the Hilbert space of solutions is 
two-dimensional, with the general solution for the zero- 



As mentioned in the Introduction, it can be expected 
on general grounds that the two-body problem with 
Coulomb interactions provides information on the more 
complicated many-body problem in graphene. This is 
especially relevant in the strong-coupling regime, as sev- 
eral works in the literature suggest the possibility of a 
new insulating phase above a certain critical coupling, 
where electron and holes would bind forming excitons 
and opening up a gap. 

It has been pointed out in the literatur o 17 ' 18 , that 
the breakdown of the Dirac vacuum in the attractive 
Coulomb impurity problem could be related to such a 
formation of excitons in graphene for strong enough cou- 
pling. We have seen in this article that this behavior is 
also present in the two-body problem, which should be 
even more relevant to the many-body physics 
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In order to understand the connection, we must notice 
that in principle the problem of an interacting electron 
and hole can be mapped into that of two attractive elec- 
trons, with similar symmetry properties. However, like 
for the Coulomb impurity problem, a more rigorous map- 
ping would be performed by considering the existence of 
the Dirac sea, which constraints, by Pauli's principle, the 
states accessible to the electron-hole pair, in analogy to 
the Cooper problem^. Such a treatment, however, would 
require to work in a different basis for which the analyt- 
ical results obtained in this paper do not hold^i, and is 
thus beyond the scope of this paper. 

We note in this regard that the breakdown of the Dirac 
vacuum in the two-body problem could be a signature of 
the excitonic instability in the many-body system. As 
we have already shown, the critical value for which the 
breakdown occurs depends on the scattering channel. For 
the most symmetric one, the s-wave, we find g c = 1, 
which should be compared to the critical values obtained 
for the Coulomb impurity problem, g^ 1 = 0.5*£~— For 
higher angular-momentum channels the critical couplings 
increase. As an example, g c — 2.24 for I = 1. How- 
ever, at low energies those higher angular momenta are 
usually less important. Hence, the s-wave critical cou- 
pling should provide an educated guess of the correspond- 
ing value for the expected many-body instability. Re- 
markably, the critical values obtained so far in the the- 
oretical literature are close to the value predicted here 
for the two-body problem. Monte Carlo calculations in 
the lattice give a critical coupling g^f ci ~ l.l l 12 ' 21 and 
gMC2 ^ 2 66, 22 depending on the model used to simu- 
late graphene electrons. Renormalization Group calcu- 
lations yield g^ G ~ 0.833. 23 Finally, a variational ap- 
proach to the excitonic condensate has been recently re- 
ported to show a transition above the critical coupling 
gvar ^ i 13,24 rp^g two-body problem with Coulomb in- 
teractions of strength above the critical coupling has not 
been addressed in this paper. As mentioned above, from 
the study of the Coulomb impurity problem it can be ex- 
pected that the instability, which produces a reorganiza- 
tion of the Dirac vacuum, leads to a non-linear screening 
of the interactions^ This effect should be analyzed care- 
fully in order to establish its connection with a possible 
formation of excitons, where it could happen that still 
the presence of the Dirac sea as a constraint is necessary 
in order to produce a bound state. After all, the results 
presented in this paper do not shed sufficient light on the 
consequences of the actual many-body instability. 

Regarding the spin degree-of-freedom, as discussed in 
section II, if both electrons belong to the same valley 
the s-wave channel would correspond to a triplet state 
in the spin sector. This fact may be highly relevant for 
the study of the excitonic instability in the presence of 
an external magnetic field. 

There is a second aspect of the two-body problem 
in graphene that could have consequences on the more 
complicated many-body problem: the influence of zero- 
energy states for I ^ angular momentum channels. As 



we have seen, in those cases where the kinetic energy van- 
ishes at some point (positive total energy and repulsive 
potential, or negative total energy and atractive poten- 
tial), zero-energy states induce singularities in the wave 
function which translate into an increasing probability 
of finding the particle at the classical turning point ro, 
when the critical coupling g c is approached. 

Let us discuss the role of the carrier density in this 
scenario. We only consider couplings below the criti- 
cal one in case of Coulomb interactions, since we ex- 
pect a description in terms of linear screening to be 
applicable. In doped graphene, electrons at the Fermi 
surface have an energy Ep — vpkp — vp^im/Ns) 1 / 2 , 
where n is the electron density in the upper cone and 
N s = 4 the valley and spin degeneracy. This defines 
a classical return distance ro = g/Ep cx n -1 / 2 for 
the Fermi surface electrons at which density correlation 
should peak. On the other hand, the static screening 
of the Coulomb interaction in doped graphene is char- 
acterized by the Thomas-Fermi (TF) screening length 
Atf = g -2 (47m./V s )~ 1 / 2 ) 2 £~— which shows a similar den- 
sity dependence, namely, Atf oc rT 1 ! 2 . Thus the ratio 
between the classical return and screening distances is 
independent of the density: to/Atf = N s g 2 . For many 
cases we expect to/Atf > 1, which places r beyond the 
screening length, i.e. where the bare Coulomb interac- 
tion, for which ro has been calculated, does not hold. 
Naively this might invalidate the physics associated to 
zero-energy states, which is expected to occur at r = ro. 
However, it is easy to note that the density correlation 
peaks have to be a robust feature of the many-body prob- 
lem. 

We have seen that zero-energy states intervene at the 
point where v(r) = E. It is quite reasonable to assume 
that, in a many-body context, that condition must be 
replaced by v SCI (r ) — Ep, which defines the classical 
return distance ro for the electron gas if v SCI (r) is the 
screened Coulomb interaction potential. Within the TF 
approximation, the screened potential has the formS^ 
^scr(fo) = (e 2 / V) F (r / Atf) , where F(x) is a monotoni- 
cally decreasing function satisfying F{x) ~ 1 for x <C 1 
and F(x) ~ x~ 2 for x ^> 1. Dimensional analysis shows 
that the dressed ro also scales like n -1 / 2 , which suggests 
that zero-energy states play a role even in the presence 
of screening. 



VIII. CONCLUSIONS 

The study of two interacting Dirac electrons in 
graphene has led us to unveil intriguing properties of 
charge carriers in this novel material. On the one hand, 
due to the chiral nature of the low-energy electrons, the 
center-of-mass and relative coordinates are coupled even 
in the presence of central potentials. This precludes a 
simple decomposition in terms of an effective one-body 
scattering problem. However, in the case of zero total 
momentum, the two-body problem can be mapped into 
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that of a single particle in the Sutherland lattice^. 

Zero-energy states turn out to play a pivotal role in the 
scattering processes, changing the matching conditions 
in the simple case of a step potential, and introducing 
singularities in the wave function for general potentials, 
including Coulomb. 

The case of Coulomb interaction is most relevant for 
the analysis of strong coupling instabilities. The rea- 
son is, electrons in weakly doped graphene have poor 
screening properties that, unlike in the conventional two- 
dimensional electron gasi22, are expected to preserve the 
long-range tail of this potential. Although the problem 
cannot be exactly mapped into the Coulomb impurity 
problem, widely studied in the literature, it still shows 
similar features such as the existence of a critical coupling 
above which wave functions become ill-defined, a likely 
signature of the Dirac vacuum breakdown. In a many- 
body context, this could signal the formation of a new 
insulating phase characterized by electron-hole pairing. 
An analysis of the s-wave scattering channel for K = 
gives a critical coupling for the instability of g c = 1, in 
rather good agreement with the critical values obtained 
in theoretical studies of the full many-body problem. We 
may also note that, due to the symmetry properties of 
the wave functions, this I = channel is accompanied by 
a spin triplet state if, as assumed throughout this paper, 
both particles belong to the same valley. 

The Coulomb potential shows other interesting proper- 
ties. We have shown that the effect of zero-energy states 
in I =/= angular channels is that of partially localizing 
the electron density near the classical turning point. The 
degree of localization increases dramatically as the crit- 
ical coupling g c is approached. Quite generally, we may 
expect the novel features found in the two-body prob- 
lem to have wide ranging implications on the many-body 
problem in graphene lattices. 
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APPENDIX A: Derivation of the matching 
conditions 

As mentioned in the main text (section IV), the so- 
lutions located at the point v(r) = E can induce dis- 
continuities in the wave function. Physically, this can 
be understood in terms of localized states that live in 
this region and which are built from the complete set of 



polynomical solutions given in (|24[) . Let us develop the 
argument in detail. 

For greater clarity, we may modify the step potential 
near the point ro where v(ro) = E in such a way that 
v(r) = E [i.e. e(r) = 0] in the vicinity of ro, namely, for 
ro — 6 < r < ro + 5, where at the end of the calculation 
S — > 0. From Eqs. (fT2"j) and flT31 in the main text, it 
follows that d r 4>2 =0 and, if I ^ 0, <f>2 = 0. On the other 
hand, Eq. (|15[) becomes, in that small interval, 

rdrfa - h) = (I - 1)^1 + {I + 1)03 • (33) 

The existence of zero-energy states [see Eq. ([M]) ] al- 
lows us to introduce functions of arbitrarily high slope 
in the small interval of length S. We adopt the simplest 
ansatze for the two components: 

Mr) = a + b(r-r ) (34) 
fo{r) = c + d(r-r ). (35) 

In the slightly modified potential, both components must 
be continuous everywhere, so we may impose 

Mro±S) = a±bS = ^ I { II (r ) (36) 
fo(r ±6) = c±dS = 4' II (r ). (37) 

As a result, if A<pi = 0f / (r o ) - $ (r ), 

6 = A0i/2<5, d = A<j> 3 /25. 

If we allow for nonzero discontinuities, A0i ^ 0, we con- 
clude that, for 6—^0, both (f>[ — b and (f)' 3 — d tend to 
infinity in magnitude. Thus, Eq. (|33p can be approxi- 
mated as 

d r Mr) = drhir) , (38) 

i.e. b = d and thus 

A<f>i = A</- 3 = -2C (39) 

The upshot is that, thanks to the existence of zero- 
energy states in the immediate vicinity of ro, a new pa- 
rameter emerges that allows for a discontinuity in the 
components </>i and <f> 3 . The parameter C is thus ad- 
justed to render the matching problem well determined. 

Interestingly, if one were to perform a similar analysis 
to the one-body problem of a step potential impurity, one 
would introduce a similar ansatz for the (only existing) 
two components of the problem. Zero-energy states could 
in principle also play a role in the vicinity of the point 
analogous to ro- However, we find that a linear ansatz 
similar to that considered above would lead to a zero 
slope. In other words, even allowing for the existence of 
zero-energy states, the wave function remains continuous 
at all points, including r . We could say that zero-energy 
states do not intervene because they are not necessary, 
and this is so because, unlike in the two-body problem, 
the matching problem is well defined from the start. 
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Once we have taken S — >■ and accepted that the 
abrupt change of sign of e(r) at r — r leads to identical 
discontinuities in <fii and 4>3 while keeping 4>2 continuous, 
we may derive, from the general relations in section III.B, 
a few more conclusions on the behavior of the solutions 
around the step. 

From the fact that e and <j>x,4>3 are bounded, it fol- 
lows from Eqs. (TT21 and (fT5|) that 4>2 and d r <p2 are also 
bounded. If we integrate Eq. (|14[) in an infintesimally 
small region around the step, we conclude 

A(e(& + &)) = , (40) 

where A means total variation across the abrupt step. 
If we combine this result with Eq. (|39|) . we conclude 
that the common discontuity of cf>x and </>3 is directly 
determined by the step discontinuity in the potential 
(Ae = —wo). Therefore Eq. (|4T))) implicitly yields the 



discontinuity C which is needed to allow <f>2 to be contin- 
uous. Specifically, we obtain 

c = \(i- l -^-)(ti + 4) , (41) 

where e 1 = E — v Q < and e 11 = E > 0. 

From Eq. (fT3")) it follows that d r (f>2 experiences a dis- 
continuity across the step which closely follows the dis- 
continuity of e(r), given that <pi — </> 3 is continuous. We 
also note from Eq. (fT2|) that, for I ^ 0, 4>2 goes quickly 
through zero as e(r) becomes at r = ro- However, it 
recovers quickly from that sharp dip to become globally 
continuous across the step, as can be inferred from (|12p 
and (|40p . By contrast, when I — 0, cf>2 remains strictly 
continuous across the step. 



I 

APPENDIX B: Analytical solution of the short-range interacting problem for zero center-of-mass momentum 



We start from the scattering problem sketched in Fig. [2j The two electron problem is written in terms of an 
effective single-electron radial equation in the case K = 0. The energy of the incident pair is E < vq. For r < ro, only 
solutions non-singular at the origin are valid, while for r > ro, a general solution is a linear combination of incoming 
and outgoing wave functions. Hence we have, up to a normalization constant [see Eqs. (j2"0)) - (l2"3")) ]. 



;Ji-i(kir) 
-Ji(kir) 
\Ji+i{kir) 



t 1 





\Ji-x(knr) 




A 


Ji{k n r) 
_ \ Ji+i{k u r) _ 


+ B 







§Ki_i(fc//r) 
Yi{k u r) 
hi(knr) 



(42) 



where the coefficients of the wave function are those of positive energy for region I and those of negative energy for 
region II. Moreover, in Eq. (|l2"]>. ki = (v - E)/2 and k u = E/2, 

As already seen in Appendix A, both solutions must be matched at r = ro with a matching condition that includes 
an arbitrary coefficient, say C, to be adjusted. The system of equations reads now: 



\Ji-i(k H r a ) ^_i(A//r ) 2 
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_ ij/ + i(fc/r ) _ 



(43) 



which can be solved by using Cramer's method. Invoking Bessel function properties, the coefficients are found to be 
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Ji(kir )-—Yi(k H ro) 
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- ^r-Yi(kiir )-^-Ji(kir ) 
ki aro 
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ki cltq 



(44) 
(45) 
(46) 



These analytical expressions reproduce the numerical results obtained by discretizing the differential equations, in- 
cluding the magnitude of the jump, — 2C, and Eq. (|4"Tj) from Appendix A. 



APPENDIX C: Analytical solution of the s-wave channel for the Coulomb interacting problem 



We start from the system of differential equations given in (|27p . The s-wave channel corresponds to I — 0. In 
this case, <j)\ — —03, and the system reduces to one of only two components, with a structure resembling that of the 
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Coulomb impurity problem. We define: 



Q\ = 4>\ - 

A 1 f. 

Ql = 01 + ^ 



which fulfill the following coupled differential equations: 



(pa p + 7 -z|)Q 1 + ^ = 

(pa p -p + 7 + l |)Q 2 + ^ i = o 



The solutions are given by Kummer functions^: 
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Qi=Ci^(7-i|,27+l;p) 
? 2 = C 2 J-(7 + l-i| ! 2 7 +l; /0 ) 



By using the property ^-"(a, b; 0) = 1 and the limit p — >• of the system of equations ([STj) . we obtain the ratio: 

C 2 



w 2 n / ■ / \ — ? arctan ^- 

C21 = 77- = - 2 (7 - *.9/ 2 ) = e 



(47) 
(48) 

(49) 

(50) 
(51) 

(52) 

(53) 
(54) 

(55) 



Hence the solution is: 



-(iEr) 



7-1/2 



e 2 



J"(7 - if ,2 7 + 1; iEr) + c 21 T(j + 1 - if ,2 7 + 1; i£V) 
2i.F( 7 - if , 2 7 + 1; i£V) - 2ic2iJ r (7 + 1 - if , 2 7 + 1; i£V) 
-^(7 - if , 2 7 + 1; iEr) - c 2i -F(7 + 1 - if , 2 7 + 1; iEr) 



(56) 



up to an overall normalization constant that can be deter mined by matching the solution to the r — > 00 limit. 
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